0.2.1 Nos packages

Voici la liste des packages que nous allons utiliser pour cet atelier :

## Pour le chargement des données, leur visualisation et manipulation
library(here)
library(dplyr)
library(DT)
library(table1)

## Pour la réalisation des tests et analyse statistiques
library(correlation)
library(parameters)
library(performance)
# Ces trois packages peuvent être appellés en une fois avec library(easystats). Easystats est une suite de packages tournée vers la réalisation de statistique et offre un environnement de travail trè sdocumenté et intéressant.
library(car)


## Représentations graphiques
library(plotly)
library(RColorBrewer)
library(ggplot2)
library(ggraph)
library(corrplot)
library(GGally)
library(gtsummary)

## Représentation cartographique
library(mapsf)

0.2.2 Nos données

0.2.3 Chargement des données sur le prix de l’immobilier par EPCI

csv_path <- here("data", "donnees_standr.csv")
immo_df <- read.csv2(csv_path)

#Pour visualiser les données dans le doc
datatable(head(immo_df, 10))

0.2.4 Chargement des données géographiques : les EPCI de France métropolitaine

Ces données proviennent de la base ADMIN-EPXRESS-COG de l’IGN, édition 2021. On les charge dans un format sf, ce qui nous permet d’utiliser le package mapsf [@mapsf] pour les prévisualiser :

library(sf)
library(mapsf)
shp_path <- here("data", "EPCI.shp")
epci_sf <- st_read(shp_path)

Reading layer `EPCI' from data source 
  `C:\Users\Lecampion\UMR_PASSAGES\AP1_ANF_FOSS\AP1_ANF_FOSS\data\EPCI.shp' 
  using driver `ESRI Shapefile'
Simple feature collection with 1242 features and 4 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: 124110 ymin: 6049642 xmax: 1242428 ymax: 7110453
Projected CRS: RGF93 v1 / Lambert-93

# pour voir les données géographiques
mf_map(x = epci_sf)

# et la table attributaire correspondante, avec le package DT
datatable(head(epci_sf, 10))

0.2.5 Jointure des données géographiques et tabulaires

Les 2 données n’ont pas le même nombre de lignes :

nrow(immo_df)

[1] 1223

nrow(epci_sf)

[1] 1242

Pour voir les EPCI qui n’ont pas de correspondance dans le tableau de données :

# l'option all.x = TRUE permet de garder toutes les lignes de epci_sf,
# même celles qui n'ont pas de correspondance dans immo_df
data_immo <- merge(x = epci_sf, y = immo_df, by.x = "CODE_SIREN", by.y = "SIREN", all.x = TRUE)
nrow(data_immo)

[1] 1242

# on peut filtrer les données de la jointure pour ne voir que les epci n'ayant pas de correspondance dans le tableau immo_df
datatable(data_immo[is.na(data_immo$prix_med),])

mf_map(x = data_immo)
mf_map(x = data_immo[is.na(data_immo$prix_med),], col = 'red', add = TRUE)

On refait la jointure en ne gardant que les EPCI ayant une correspondance dans le tableau de données (on choisit donc de laisser de côté les EPCI non listés dans le fichier CSV)

data_immo <- merge(x = epci_sf, y = immo_df, by.x = "CODE_SIREN", by.y = "SIREN")
nrow(data_immo)

[1] 1223

datatable(head(data_immo, 10))

0.4.1 Micro rappel théorique

La corrélation permet d’étudier le lien, la relation entre deux variables.

La corrélation repose sur la covariance entre les variables. Quand 2 variables covarient, un écart à la moyenne d’une variable est accompagné par un écart dans le même sens (corrélation positive) ou dans le sens opposé (corrélation négative) pour le même individu de l’autre variable. Dit autrement, lorsque deux variables covarient, pour chaque valeur qui s’écarte de la moyenne, on s’attend à trouver un écart à la moyenne pour l’autre variable.

Une corrélation est nécessairement comprise entre 1 et -1.

Un coefficient de 1 indique une corrélation positive parfaite entre les deux variables. Quand la variable \(x\) augmente la variable \(y\) augmente dans la même proportion. À l’inverse, un coefficient de – 1 indique une corrélation négative parfaite : lorsque la variable \(x\) augmente, la variable \(y\) diminue dans la même proportion. Dans les deux cas, les points tombent parfaitement sur la droite. Un coefficient de 0 indique qu’il n’y a aucune relation entre les deux variables. Ainsi, la variation de l’une n’est aucunement associée à la variation de l’autre.

Par exemple, ici le % de petits logements est corrélé positivement au taux de logements suroccupés : si l’une de ces 2 variables augmente, l’autre aussi. A l’inverse, le % de maisons est corrélé négativement au taux de logement suroccupés : si l’une de ces 2 variables augmente, l’autre diminue.

Pour évaluer l’importance d’une corrélation on utilise souvent la convention suivante :

Interprétation d’une corrélation

L’analyse de corrélation présente plusieurs intérêts. Le premier étant bien sur d’étudier les relations entre les variables, ce qui en soit peut s’avérer très intéressant. Ensuite, la conception d’un modèle doit absolument être le fruit d’une réflexion portant sur le choix des variables indépendantes (explicatives) et le choix de la méthode de régression. Et pour définir un modèle de régression certaine règles doivent (dans un monde parfait) être respectées. L’étude des corrélations peut donc apporter une aide précieuse dans cette réflexion. Elle pourra nous aider dans le choix des variables à intégrer au modèle, et dans le même temps à vérifier certaines des conditions de réalisation de notre régression.

Ainsi, une analyse des corrélations pourra vérifier :

• l’existence d’un lien entre les variables indépendantes (VI) et la variable à étudier (VD). En effet dans une régression linéaire, il est nécessaire d’avoir une relation linéaire entre la VD et les différentes VI.
• la multicolinéarité des variables indépendantes. Les corrélations ne doivent pas être trop fortes entre nos VI. Un coefficient > 0.7 en valeurs absolues doit entraîner la suppression des variables concernés. Cela peut aussi être vérifié très efficacement avec le VIF (Variance Inflation Factor) mais peut se faire seulement après avoir lancé notre modèle.
• L’absence de corrélation entre les variables explicatives du modèle et les variables externes. En effet, les variables d’influence doivent être incluses dans le modèle (sauf dans le cas où cela induirait une trop grande multicolinéarité).

Pour calculer nos corrélation nous utiliserons le package bien nommé correlation. Ce package présente de nombreux avantages, comme plusieurs manière de présenter nos résultats (base de données, matrice…). Il intégre par défault des tests de corrections qui permettent des calculs plus robustes, il est compatible avec la méthode bayesienne, propose un grand nombre de méthode différentes…

0.4.2 Réalisation des correlations avec R

L’analyse se lance simplement avec ces lignes de codes:

library(correlation)

data_cor <- immo_df %>% select(-SIREN)
immo_cor <- correlation(data_cor)

Un des avantages de l’analyse de corrélations c’est que les résultats sont facilement interprétables et accessible de beaucoup de manière différentes

# En mode base de données
datatable(head(immo_cor, 10))

# En matrice des corrélations
summary(immo_cor)

# Correlation Matrix (pearson-method)

Parameter           | part_cadre_profintellec_nbemploi | part_agri_nb_emploi | part_log_suroccup | med_niveau_vis | dens_pop | perc_tiny_log | perc_maison | perc_log_vac
-------------------------------------------------------------------------------------------------------------------------------------------------------------------------
prix_med            |                          0.55*** |            -0.45*** |           0.60*** |        0.59*** |  0.52*** |       0.49*** |    -0.60*** |     -0.68***
perc_log_vac        |                         -0.31*** |             0.39*** |          -0.28*** |       -0.46*** | -0.21*** |      -0.19*** |     0.35*** |             
perc_maison         |                         -0.58*** |             0.54*** |          -0.79*** |       -0.24*** | -0.46*** |      -0.75*** |             |             
perc_tiny_log       |                          0.75*** |            -0.51*** |           0.87*** |        0.22*** |  0.62*** |               |             |             
dens_pop            |                          0.59*** |            -0.29*** |           0.62*** |        0.18*** |          |               |             |             
med_niveau_vis      |                          0.43*** |            -0.33*** |           0.15*** |                |          |               |             |             
part_log_suroccup   |                          0.65*** |            -0.43*** |                   |                |          |               |             |             
part_agri_nb_emploi |                         -0.57*** |                     |                   |                |          |               |             |             

p-value adjustment method: Holm (1979)

On peut aussi visualiser les corrélations de nombreuses manières différentes, il y a le classique scatter plot pour visualier la relation entre deux variables mais pour représenter notre matrice la figure classique est le corrélogramme :

1. Méthode simple avec le package see

library(see)
immo_cor %>%
  summary(redundant = FALSE) %>%
  plot(type="tile", show_labels =TRUE, show_p = TRUE, digits = 1, size_text=3) +
  see::theme_modern(axis.text.angle = 45)

2. Méthode plus complexe avec le package corrplot mais permet d’avantage de personnalisation

#création d'une matrice redondante
mat_cor_comp <- summary(immo_cor, redundant = TRUE)
# Nom des lignes = valeurs de la première colonne ("Parameter")
rownames(mat_cor_comp ) <- mat_cor_comp[,1]
# Transformation du data.frame en objet matrice (+ suppression première colonne)
mat_cor<- as.matrix(mat_cor_comp[,-1])
# Calcul du nombre total d'individus
nb <- nrow(data_cor)
# Calcul des matrices de p-values et des intervalles de confiance
p.value <- cor_to_p(mat_cor, n = nb, method = "auto")
# Extraction de la matrice des p-value uniquement
p_mat <- p.value$p

library(corrplot)
library(RColorBrewer)

corrplot(mat_cor, 
         p.mat = p_mat, 
         type = "upper", 
         order = "hclust", 
         addCoef.col = "white", 
         tl.col = "gray",
         number.cex = 0.5,
         tl.cex= 1,
         tl.srt = 45, 
         col=brewer.pal(n = 8, name = "PRGn"), 
         sig.level = 0.000001, 
         insig = "blank", 
         diag = FALSE, )

La représentation sous forme de réseau peut être très intéressante notamment quand il y a un grand nombre de variables, et donc de corrélations.

library(ggraph) # needs to be loaded

immo_cor$val <-if_else(
  immo_cor$r>0, "positif", "négatif"
)

ggraph(immo_cor, layout = "stress") +
 geom_node_point(size = 8, color="grey")+
  geom_node_text(aes(label = name),family = "serif",repel = TRUE,  max.overlaps = getOption("ggrepel.max.overlaps", default = Inf))+
  geom_edge_link2(aes(filter = p<0.05, edge_colour= val, edge_width = abs(r), label = round(r, 2))) + # on peut indiquer sur le lien avec l'argument  label = r
  scale_edge_width_continuous(range = c(0.1,3))+
scale_edge_color_manual(
    values = c("positif" = "#1B9E77",
               "négatif" = "#D95F02")
  )

Pour vérifier la relation entre notre VD et une VI on peut classiquement faire un scatter plot, ici par exemple entre le prix médian (VD) et le pourcentage de logements vacants (VI).

# Relations bivariées - formes fonctionnelles
ggplot(data_immo, aes(x=log(perc_log_vac), y=log(prix_med))) + 
  geom_point() + geom_smooth()

Le nuage de point est l’autre représentation standard des corrélation. Il s’agit ici d’analyse la forme du nuage pour avoir une information sur la corrélation. Cette relation n’est pas la simple à lire mais ce qu’il faut retenir c’est plus vos points sont regroupé autour de la droute de régession plus il y a de la chance que la corrélation soit significative. La forme du nuage nous donnera le signe de la corrélation. Un nuage qui part d’en bas à gauche vers en haut à droite aura une corrélation positive et un nuage qui va d’en haut à gauche vers en bas à droite aura une corrélation négative.

Mais ce qui serait plus intéressant serait une représentation de notre VDvis à vis de toutes nos VI. On peut utiliser alors un diagramme en barres.

library(corrr)

# Matrice de corrélation -> data.frame de corrélation
hist_cor <- as_cordf(mat_cor)

# Sélection des corrélations de la variable "year"
hist_cor <- focus(hist_cor, prix_med)
# Si la valeur est positive  = TRUE, sinon FALSE
hist_cor$correlation <- hist_cor$prix_med > 0

ggplot(hist_cor, aes(x = reorder(term, prix_med), y = prix_med, fill = correlation)) +
      geom_bar(stat = "identity", colour = "black") +
      ylab("Corrélation") +
      xlab("Variables") +
      scale_fill_discrete(labels = c("Négative", "Positive")) +
      theme_modern(axis.text.angle = 45) 

0.5.1 Micro rappel sur la régression

La régression linéaire est un des modèles les plus utilisés en SHS, elle peut être simple (une seule variable explicative) ou multiple (plusieurs VI). Le principe n’est en réalité pas si complexe. La régression linéaire consiste à modéliser la covariation entre une variable à expliquer et une ou des variables explicatives.

Pour ce faire, le modèle de régression va chercher à estimer les termes de l’équation de la droite de régression entre la variable à expliquer et la variable explicative. Cette équation va prendre la forme \(f(x)= ax + b + e_i\). Équation qui ressemble à la fonction affine qui a dû hanter vos nuits de collégien.

Si on reprend l’exemple précédent où l’on cherche à modéliser la covariation entre le prix médian et le pourcentage de logement vacant, l’équation de notre droite de régression ressemblerait à : \(prixmed_i = a*perclogvac_i + b + e_i\).

Où \(a\) est le coefficient associé à la variable du pourcentage de logements vacants, \(b\) est la constante qui apparaîtra sous le nom d’intercept dans les résultats et enfin \(e\) qui lui correspond aux résidus. Les résidus étant ce qui incarne l’écart au modèle.

Traditionnellement on va faire usage d’une régression linéaire lorsque l’on veut prédire les valeurs de notre VD dans les cas où elle n’aurait pas été mesurée, ou bien si l’on souhaite comprendre les relations statistiques entre les variables. Comment le sprédicteurs vont expliquer notre variable dépendante.

Pour lancer notre modèle de régression linéaire, on va utiliser la fonction lm() dont les 2 lettres sont l’acronyme de linear model :

Sur R la réalisation d’un modèle de régression linéaire se fait très simplement :

# Dans le fonctionnement sur R il est important de stocker la régression dans un objet.
mod.lm <- lm(formula = prix_med ~ perc_log_vac + perc_maison + perc_tiny_log + dens_pop + med_niveau_vis + part_log_suroccup + part_agri_nb_emploi + part_cadre_profintellec_nbemploi, 
             data = data_immo)

# On affiche les principaux résultats avec la fonction summary
summary(mod.lm)

Call:
lm(formula = prix_med ~ perc_log_vac + perc_maison + perc_tiny_log + 
    dens_pop + med_niveau_vis + part_log_suroccup + part_agri_nb_emploi + 
    part_cadre_profintellec_nbemploi, data = data_immo)

Residuals:
    Min      1Q  Median      3Q     Max 
-1566.8  -220.2   -27.2   174.0  3277.3 

Coefficients:
                                 Estimate Std. Error t value Pr(>|t|)    
(Intercept)                      1592.873     11.204 142.171  < 2e-16 ***
perc_log_vac                     -287.337     14.134 -20.330  < 2e-16 ***
perc_maison                      -128.255     20.041  -6.400 2.22e-10 ***
perc_tiny_log                    -261.556     27.934  -9.363  < 2e-16 ***
dens_pop                          173.942     15.272  11.389  < 2e-16 ***
med_niveau_vis                    288.017     13.860  20.781  < 2e-16 ***
part_log_suroccup                 388.982     27.352  14.221  < 2e-16 ***
part_agri_nb_emploi               -20.785     15.059  -1.380    0.168    
part_cadre_profintellec_nbemploi   -7.841     19.904  -0.394    0.694    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 391.8 on 1214 degrees of freedom
Multiple R-squared:  0.7784,    Adjusted R-squared:  0.7769 
F-statistic: 532.9 on 8 and 1214 DF,  p-value: < 2.2e-16

Une des fonction les plus utile dans tout l’écosystème R est la fonction summary(), ici pour notre régression elle va nous fournir tous les élément utile à l’interprétation. Voci donc ce qu’elle nous sort :

• Call : la fomule du modèle ;
• Residuals : des statistiques descriptives des résidus ;
• Coefficients : un tableau à deux entrées où les lignes correspondent aux coefficients associés aux variables explicatives, et les colonnes, dans l’ordre, à l’estimation du coefficient, l’écart-type estimé, la valeur du test de Student de nullité statistique du coefficient et enfin la p-value associé à ce test, suivie d’un symbole pour lire rapidement la significativité;
• Signif. codes : les significations des symboles de niveau de significativité ;
• Residual standard error : estimation de l’écart-type de l’aléa et degré de liberté ;
• Multiple R-squared : coefficient de détermination ;
• Adjusted R-squared : coefficient de détermination ajusté ;
• F-statistic : valeur de la statistique de Fisher du test de significativité globale, ainsi que les degrés de liberté et la p-value associée au test.

Pour visualiser les résultats des tests statistiques il existent un grand nombre de package et de fonctions. On peut citer par exemple la fonction tbl_regression du package gtsummary ou la fonction model_parameters du package parameters.

library(gtsummary)
mod.lm %>%
  tbl_regression(intercept = TRUE)

library(parameters)
model_parameters(mod.lm) |> print_md()

On peut visualiser graphiquement les coefficients des variables explicatives, ce qui peut s’avérer très utile :

GGally::ggcoef_model(mod.lm)

Faire une régression linéaire sur R c’est finalement très simple ! Mais le vrai enjeux c’est l’interprétation des résultats !

0.5.2 Interpréter les résultats

Pour interpréter les résultats plusieurs éléments fournis par la commande summary(mod.lm) sont importants.

Résultats de la commande summary(mod.lm)

Commençons par les informations qui concernent notre modèle dans sa globalité. D’abord l’information fournie par Adjusted R-squared, il s’agit du \(R^2\) qui est le coefficient de détermination. Il est ici de 0.77 ce qui peut se traduire par 77% de la variation du prix médian est expliqué par notre modèle. Juste en-dessous, on a le p-value associé à notre modèle. S’il est inférieur à \(0.05\) on peut considérer que le modèle est statistiquement significatif.

Ensuite, regardons ce qu’il se passe au niveau des variables explicatives :

• La première colonne appelée estimate est le coefficient de régression associé à la variable explicative. Le signe va être très important car il va donner la direction de la relation (exactement comme pour les corrélations). Ici pour le pourcentage de logements vacants il est de \(-287\) ce qui veut dire que lorsque ce pourcentage augmente de \(1\%\) alors le prix médian baisse de \(287€\)

• La colonne \(Pr(>|t|)\) correspond à la p-value associé à ce résultat. Une fois encore si elle est inférieure à \(0.05\) alors on peut le considérer comme statistiquement significatif. Dans notre cas on peut dire que la part d’agriculteurs, de cadres et professions intellectuelles dans le nombre d’emploi n’ont pas un effet significatif sur le prix médian du logement.

Finalement tout va pour le mieux dans le meilleur des mondes ! On a donc pu réaliser notre régression et interpréter les résultats. Mais en réalité, nous avons pas abordé le point le plus sensible lorsque l’on réalise un modèle de régression linéaire. Notre modèle est-il viable ? Les conditions nécessaire d’application sont-elle respecté ? Et c’est bien là le plus dure dans la régression.

0.6.1 La Multicolinéarité

Un des enjeux les plus important dans le cadre de régression multiple est de vérifier la multicolinéarité entre les variables explicatives. Le risque d’une trop grande colinéarité est de biaiser le modèle et notamment de biaiser les estimations des erreurs type des coefficient de régression (et donc les t-value et p-value).

L’analyse préalable des corrélations donne déjà de forte présemption sur un risuqe de multicolinéarité mais attention il ne faut pas confondre la multicolinéarité avec une “forte corrélation brutes” entre les prédicteurs. La multicolinéarité entre deux prédicteurs peut survenir lorsqu’une troisième variable non observée a un effet causal sur chacun des deux prédicteurs associés au résultat. Dans de tels cas, la relation réelle qui compte serait l’association entre la variable non observée et le résultat.

La VIF (Variance Inflation Factor) est une très bonne méthode pour vérifier et peut être confirmer les risques de multicolinéarité. Elle suppose simplement d’avoir estimé un premier modèle pour être utilisée. Cela se fait très simplement sur R, je vus montre plusieurs manière de le faire à vous de choisir celle que vous préférez.

# Avec la librairie "car"
library(car)

vif(mod.lm)

                    perc_log_vac                      perc_maison 
                        1.577984                         3.204058 
                   perc_tiny_log                         dens_pop 
                        6.221631                         1.864236 
                  med_niveau_vis                part_log_suroccup 
                        1.532403                         5.977951 
             part_agri_nb_emploi part_cadre_profintellec_nbemploi 
                        1.812419                         3.162576 

# Avec la librairie "performance" de la suite easystats qui propose beaucoup de choses globalement pour vous faciliter la vie si vous faites de la modélisation statistique. Toutefois attention ce sont plutôt des psychologues qui en sont à l'origine et ce sont donc plutôt leurs normes disciplinaires.

library(performance)

check_collinearity(mod.lm)

# Check for Multicollinearity

Low Correlation

                             Term  VIF   VIF 95% CI Increased SE Tolerance
                     perc_log_vac 1.58 [1.47, 1.71]         1.26      0.63
                      perc_maison 3.20 [2.92, 3.53]         1.79      0.31
                         dens_pop 1.86 [1.73, 2.03]         1.37      0.54
                   med_niveau_vis 1.53 [1.43, 1.66]         1.24      0.65
              part_agri_nb_emploi 1.81 [1.68, 1.97]         1.35      0.55
 part_cadre_profintellec_nbemploi 3.16 [2.89, 3.48]         1.78      0.32
 Tolerance 95% CI
     [0.58, 0.68]
     [0.28, 0.34]
     [0.49, 0.58]
     [0.60, 0.70]
     [0.51, 0.60]
     [0.29, 0.35]

Moderate Correlation

              Term  VIF   VIF 95% CI Increased SE Tolerance Tolerance 95% CI
     perc_tiny_log 6.22 [5.62, 6.90]         2.49      0.16     [0.14, 0.18]
 part_log_suroccup 5.98 [5.41, 6.62]         2.44      0.17     [0.15, 0.19]

# On peut aussi directement l'ajouter au résumé des coefficient obtenu avec gtsummary

library(gtsummary)

mod.lm %>%
  tbl_regression(intercept = TRUE) %>% add_vif()

On peut facilement représenter graphiquement les scores de VIF.

# Avec la librairie "car"
library(car)

score_vif <- vif(mod.lm)

barplot(score_vif, main = "VIF Values", horiz = TRUE, col = "steelblue", las=2)
#ajout du seuil de 4
abline(v = 4, lwd = 3, lty = 2)
# et de la limite de 3
abline(v = 3, lwd = 3, lty = 2)

# Avec la librairie performance, mais attention ici pour le VIF il 'agit d'un seuil qui est représenté à voir s'il est compatible avec la pratique de votre discipline

colinearite <-  check_collinearity(mod.lm)
plot(colinearite)

Comment gérer la multicolinéarité ?

Nous avons en effet une forte multicolinéarité entre certaines de nos variables explicatives. Selon le VIF nous devrions donc relancer notre modèle sans les variables fortement colinéaires, c’est à dire sans le pourcentage de logements sur-occupés et sans le pourcentage de petits logements.

Normalement, il est d’usage de retirer ces variables de notre modèle. Voyons ce qui se passe. Pour commencer, on peut retirer du modèle la variable ayant le VIF le plus élevé à savoir le pourcentage de petits logements.

mod.lm2 <- lm(formula = prix_med ~ perc_log_vac + perc_maison + dens_pop + 
    med_niveau_vis + part_log_suroccup + part_agri_nb_emploi + 
    part_cadre_profintellec_nbemploi, data = data_immo)

summary(mod.lm2)

Call:
lm(formula = prix_med ~ perc_log_vac + perc_maison + dens_pop + 
    med_niveau_vis + part_log_suroccup + part_agri_nb_emploi + 
    part_cadre_profintellec_nbemploi, data = data_immo)

Residuals:
    Min      1Q  Median      3Q     Max 
-1502.2  -226.3   -42.3   173.5  3535.9 

Coefficients:
                                 Estimate Std. Error t value Pr(>|t|)    
(Intercept)                      1592.549     11.597 137.329  < 2e-16 ***
perc_log_vac                     -328.301     13.911 -23.601  < 2e-16 ***
perc_maison                      -100.010     20.507  -4.877 1.22e-06 ***
dens_pop                          161.783     15.750  10.272  < 2e-16 ***
med_niveau_vis                    280.586     14.322  19.591  < 2e-16 ***
part_log_suroccup                 237.069     22.792  10.401  < 2e-16 ***
part_agri_nb_emploi                 2.696     15.370   0.175    0.861    
part_cadre_profintellec_nbemploi  -77.744     19.098  -4.071 4.99e-05 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 405.5 on 1215 degrees of freedom
Multiple R-squared:  0.7623,    Adjusted R-squared:  0.761 
F-statistic: 556.8 on 7 and 1215 DF,  p-value: < 2.2e-16

vif(mod.lm2)

                    perc_log_vac                      perc_maison 
                        1.426783                         3.131469 
                        dens_pop                   med_niveau_vis 
                        1.850756                         1.527378 
               part_log_suroccup              part_agri_nb_emploi 
                        3.874580                         1.762159 
part_cadre_profintellec_nbemploi 
                        2.717634 

library(gtsummary)
mod.lm2 %>%
  tbl_regression(intercept = TRUE) %>% add_vif()

GGally::ggcoef_model(mod.lm2)

On note qu’au niveau global il y a peu de changements : le \(R^2\) a très légèrement baissé, on passe d’une variation expliquée à 77% à un taux d’explication de 76% et le modèle est toujours significatif. Les changement les plus importants se situent au niveau des effets partiels. L’effet de la part de cadre et de professions intellectuelles dans le nombre d’emploi sur le prix médian du logement est devenu significatif et on constate même que le VIF du pourcentage de logement sur-occupés est passé sous le seuil critique.

step(mod.lm2, direction = "backward")

Start:  AIC=14696.68
prix_med ~ perc_log_vac + perc_maison + dens_pop + med_niveau_vis + 
    part_log_suroccup + part_agri_nb_emploi + part_cadre_profintellec_nbemploi

                                   Df Sum of Sq       RSS   AIC
- part_agri_nb_emploi               1      5060 199817567 14695
<none>                                          199812507 14697
- part_cadre_profintellec_nbemploi  1   2725360 202537867 14711
- perc_maison                       1   3911245 203723752 14718
- dens_pop                          1  17351111 217163618 14796
- part_log_suroccup                 1  17792190 217604698 14799
- med_niveau_vis                    1  63121115 262933622 15030
- perc_log_vac                      1  91601521 291414028 15156

Step:  AIC=14694.71
prix_med ~ perc_log_vac + perc_maison + dens_pop + med_niveau_vis + 
    part_log_suroccup + part_cadre_profintellec_nbemploi

                                   Df Sum of Sq       RSS   AIC
<none>                                          199817567 14695
- part_cadre_profintellec_nbemploi  1   3211704 203029272 14712
- perc_maison                       1   4155413 203972980 14718
- dens_pop                          1  17567092 217384659 14796
- part_log_suroccup                 1  18007058 217824626 14798
- med_niveau_vis                    1  63117884 262935451 15028
- perc_log_vac                      1  94962190 294779757 15168

Call:
lm(formula = prix_med ~ perc_log_vac + perc_maison + dens_pop + 
    med_niveau_vis + part_log_suroccup + part_cadre_profintellec_nbemploi, 
    data = data_immo)

Coefficients:
                     (Intercept)                      perc_log_vac  
                         1592.55                           -327.82  
                     perc_maison                          dens_pop  
                          -99.01                            162.05  
                  med_niveau_vis                 part_log_suroccup  
                          280.55                            237.44  
part_cadre_profintellec_nbemploi  
                          -78.93  

mod.lm3 <- lm(formula = prix_med ~ perc_log_vac + perc_maison + dens_pop + med_niveau_vis + part_log_suroccup + part_cadre_profintellec_nbemploi, data = data_immo)

library(performance)

perf <- compare_performance(mod.lm, mod.lm2, mod.lm3)
print_md(perf)

library(see)
plot(perf)

0.6.2 Analyser les résidus

L’analyse des résidus est très importante car les conditions de validité d’un modèle linéaire au delà des résultats repose grandement sur les résidus. Ils permettent en outre aussi d’identifier les individus extrêmes (ou outliers).

Les 3 conditions qui concernent les résidus sont :

• Ils doivent suivre une loi normale.
• Ils ne doivent pas varier en fonction des variables explicatives. C’est l’hypothèse d’homoscédasticité, ils ont une variance homogène.
• Ils ne doivent pas être autocorrélés.

Pour obtenir les résidus :

res_modlm <- mod.lm$residuals
datatable(as.data.frame(res_modlm))

Brutes de cette manière cela ne dit pas grand chose, il est nécessaire de les visualiser. De manière générale voici les graphiques que l’on voit fréquemment concernant les résidus en R base.

par(mfrow=c(1,3))
qqPlot(mod.lm) # diagramme quantile-quantile qui permet de vérifier l'ajustement d'une distribution à un modèle théorique, ici loi normale

[1]  36 266

hist(rstudent(mod.lm), breaks = 50, col="darkblue", border="white", main="Analyse visuelle des résidus") # Histogramme pour donner une autre indication sur la normalité
plot(rstudent(mod.lm)) # un graphique pour visualiser l'homoscédasticité des résidus

Lorsque graphiquement, il est compliqué de se prononcer, il est possible d’utiliser des tests statistqique qui ont pour objet soit la vérification de la normalité ou de tester l’homoscédasticité. Ils ont cela de particulier qu’ici nous cherchons à accepter H0 et donc pour valider la normalité ou l’homoscédasticité il faut que \(p-value > 0.05\)

# Pour étudier la normalité on peut utiliser le test de Shapiro-Wilk
shapiro.test(mod.lm$residuals)

    Shapiro-Wilk normality test

data:  mod.lm$residuals
W = 0.90792, p-value < 2.2e-16

# Pour évaluer l'homoscédasticité on peut utiliser le test de Breusch-Pagan. Le package car propose une fonction pour le réaliser
ncvTest(mod.lm)

Non-constant Variance Score Test 
Variance formula: ~ fitted.values 
Chisquare = 1151.844, Df = 1, p = < 2.22e-16

Dans les deux cas il nous faut rejeter H0, les résidus n’ont donc pas une distribution normale et il y a hétéroscédasticité de la variance des résidus.

Le modèle ne peut donc pas être analysé en l’état. Le problème de l’hétéroscédasticité des résidus indique souvent un problème de spécification du modèle (par exemple une variable manquante).

Le package performance peut une fois de plus être une alternative très intéressante et nous macher le travail. On pourra aisément tester ces conditions et les visualiser de manière moins aride qu’en R base. Ces fonctions fournissent des indicateurs d’aide à la décision qui sont plus claires.

# Pour la normalité des résidus
normality <- check_normality(mod.lm)
normality

Warning: Non-normality of residuals detected (p < .001).

plot(normality)

# Pour l'homoscédasticité

heteroscedasticity <- check_heteroscedasticity(mod.lm)
heteroscedasticity

Warning: Heteroscedasticity (non-constant error variance) detected (p < .001).

plot(heteroscedasticity)

ind_ext <- check_outliers(mod.lm)
ind_ext

1 outlier detected: case 266.
- Based on the following method and threshold: cook (0.927).
- For variable: (Whole model).

plot(ind_ext, type = "dots")

0.6.3 Apparté sur les outliers

Lors de la modélisation statistique on passe énormément de temps à réfléchir aux variables, à comment ces dernières peuvent venir perturber notre modèle. Cependant, certains de nos individus peuvent eux aussi venir grandement influencer nos analyses : il s’agit des individus extrêmes ou outliers. LLes analyses que nous venons de faire nous indique globalement deux individus particulier, ici ceux ayant pour identifiant 266 et 36. Il peut dans certain cas être intéressant de supprimer ces individus et voir comment réagit le modèle.

# Pour visualiser les individus concernés
data_immo[c(36, 266),]

Simple feature collection with 2 features and 13 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: 637302.2 ymin: 6532738 xmax: 1013257 ymax: 6879128
Projected CRS: RGF93 v1 / Lambert-93
    CODE_SIREN                       ID                                    NOM
36   200023372 EPCI____0000002150235984 CC de la Vallée de Chamonix-Mont-Blanc
266  200054781 EPCI____0000002150236221               Métropole du Grand Paris
                    NATURE prix_med perc_log_vac perc_maison perc_tiny_log
36  Communauté de communes     6400   -2.3174835   -2.558919   -0.07390879
266              Métropole    10920   -0.7044846   -3.732690    6.43751449
      dens_pop med_niveau_vis part_log_suroccup part_agri_nb_emploi
36  -0.2437935      0.7430146          1.873292           -1.027161
266 22.7709258      0.9077121          7.345980           -1.068415
    part_cadre_profintellec_nbemploi                       geometry
36                        0.07587777 MULTIPOLYGON (((997329.6 65...
266                       6.30681026 MULTIPOLYGON (((666798.7 68...

# Pour relancer un nouveau modèle sans l'individu le plus extrême. Notez que l'on peut en supprimer plusieurs d'un coup avec subset =-c(36,266)
mod.lmx <- update(mod.lm, subset=-266)

# Etudier le nouveau modèle

summary(mod.lmx)

Call:
lm(formula = prix_med ~ perc_log_vac + perc_maison + perc_tiny_log + 
    dens_pop + med_niveau_vis + part_log_suroccup + part_agri_nb_emploi + 
    part_cadre_profintellec_nbemploi, data = data_immo, subset = -266)

Residuals:
    Min      1Q  Median      3Q     Max 
-1622.6  -215.0   -37.5   170.0  2767.7 

Coefficients:
                                 Estimate Std. Error t value Pr(>|t|)    
(Intercept)                      1587.246     10.526 150.797  < 2e-16 ***
perc_log_vac                     -302.182     13.317 -22.692  < 2e-16 ***
perc_maison                      -132.159     18.814  -7.024 3.57e-12 ***
perc_tiny_log                    -227.353     26.356  -8.626  < 2e-16 ***
dens_pop                           -4.725     19.978  -0.237    0.813    
med_niveau_vis                    284.947     13.012  21.899  < 2e-16 ***
part_log_suroccup                 404.159     25.701  15.725  < 2e-16 ***
part_agri_nb_emploi               -17.545     14.138  -1.241    0.215    
part_cadre_profintellec_nbemploi   23.384     18.841   1.241    0.215    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 367.8 on 1213 degrees of freedom
Multiple R-squared:  0.7823,    Adjusted R-squared:  0.7809 
F-statistic: 544.9 on 8 and 1213 DF,  p-value: < 2.2e-16

vif(mod.lmx)

                    perc_log_vac                      perc_maison 
                        1.589305                         3.168447 
                   perc_tiny_log                         dens_pop 
                        6.073006                         2.089258 
                  med_niveau_vis                part_log_suroccup 
                        1.531885                         5.726830 
             part_agri_nb_emploi part_cadre_profintellec_nbemploi 
                        1.811307                         3.111740 

# Il est possible de comparer les deux modèles et les coefficients

car::compareCoefs(mod.lm, mod.lmx, pvals = TRUE)

Calls:
1: lm(formula = prix_med ~ perc_log_vac + perc_maison + perc_tiny_log + 
  dens_pop + med_niveau_vis + part_log_suroccup + part_agri_nb_emploi + 
  part_cadre_profintellec_nbemploi, data = data_immo)
2: lm(formula = prix_med ~ perc_log_vac + perc_maison + perc_tiny_log + 
  dens_pop + med_niveau_vis + part_log_suroccup + part_agri_nb_emploi + 
  part_cadre_profintellec_nbemploi, data = data_immo, subset = -266)

                                 Model 1 Model 2
(Intercept)                       1592.9  1587.2
SE                                  11.2    10.5
Pr(>|z|)                         < 2e-16 < 2e-16
                                                
perc_log_vac                      -287.3  -302.2
SE                                  14.1    13.3
Pr(>|z|)                         < 2e-16 < 2e-16
                                                
perc_maison                       -128.3  -132.2
SE                                  20.0    18.8
Pr(>|z|)                         1.6e-10 2.1e-12
                                                
perc_tiny_log                     -261.6  -227.4
SE                                  27.9    26.4
Pr(>|z|)                         < 2e-16 < 2e-16
                                                
dens_pop                          173.94   -4.72
SE                                 15.27   19.98
Pr(>|z|)                         < 2e-16    0.81
                                                
med_niveau_vis                     288.0   284.9
SE                                  13.9    13.0
Pr(>|z|)                         < 2e-16 < 2e-16
                                                
part_log_suroccup                  389.0   404.2
SE                                  27.4    25.7
Pr(>|z|)                         < 2e-16 < 2e-16
                                                
part_agri_nb_emploi                -20.8   -17.5
SE                                  15.1    14.1
Pr(>|z|)                            0.17    0.21
                                                
part_cadre_profintellec_nbemploi   -7.84   23.38
SE                                 19.90   18.84
Pr(>|z|)                            0.69    0.21
                                                

# Si on r

perf2 <- compare_performance(mod.lm, mod.lmx)
print(perf2)

# Comparison of Model Performance Indices

Name    | Model |   AIC (weights) |  AICc (weights) |   BIC (weights) |    R2 | R2 (adj.) |    RMSE |   Sigma
-------------------------------------------------------------------------------------------------------------
mod.lm  |    lm | 18086.1 (<.001) | 18086.3 (<.001) | 18137.2 (<.001) | 0.778 |     0.777 | 390.353 | 391.797
mod.lmx |    lm | 17916.6 (>.999) | 17916.8 (>.999) | 17967.7 (>.999) | 0.782 |     0.781 | 366.405 | 367.762

plot(perf2)

Dans certains cas les différences sont mineures, ici la différence est importante, en effet, on voit qu’une VI a perdu sa significativité.

Quoi qu’il en soit c’est une opération qui doit être effectuée avec prudence, la suppression des individus pose toujours question notamment en terme de justification théorique. Il faut le faire uniquement si l’analyse des individus indique un problème important (valeurs aberrantes, inversion…)

0.6.4 L’autocorrélation des résidus

C’est la condition qui pose en généralle le plus problème, car elle implique beaucoup de choses finalement sur nos données ou notre cadre de réflexions.

Pour vérifier facilement s’il y a un quelconque risque d’auto-correlation :

check_autocorrelation(mod.lm)

Warning: Autocorrelated residuals detected (p < .001).

Cependant, de cette manière on ne sait pas du tout de quel ordre est cette autocorrélation. Nos données étant spatialisé, il serait pertinent de vérifier spécifiquement le risque d’autocorrélation spatiale.

Heureusement la géographie s’est doté d’outils pour mesurer notamment l’autocorrélation spatiale. En réalité ici nous espérons très fortement qu’il y ait une autocorrélation spatiale. Cela rendrait notre modèle linéaire classique caduc mais nous permettrait de justifier l’utilisation de la régression spatiale et donc cette seamaine d’ANF

Les deux outils connus au moins de nom par tous les géographes sont le test de Moran et celui de Geary.

Dans la littérature le test de Moran semble être préféré à celui de Geary en raison d’une stabilité plus grande, mais dans les faits les deux sont souvent associé

# test de moran des residus de la régression H0: pas d'autocorrélation spatiale
lm.morantest(model = mod.lm, 
            listw = neighbours_epci_w)

    Global Moran I for regression residuals

data:  
model: lm(formula = prix_med ~ perc_log_vac + perc_maison +
perc_tiny_log + dens_pop + med_niveau_vis + part_log_suroccup +
part_agri_nb_emploi + part_cadre_profintellec_nbemploi, data =
data_immo)
weights: neighbours_epci_w

Moran I statistic standard deviate = 21, p-value < 2.2e-16
alternative hypothesis: greater
sample estimates:
Observed Moran I      Expectation         Variance 
     0.359707408     -0.003517610      0.000299161 

# Test de Geary H0 pas d'autocorrélation.
#  Attention : Pour avoir le  coefficient il faut faire 1-"Résultat test de Geary" (soit ici le coefficient est 0.67)
# Le coefficient de Geary s'étend de 0 à 2, 1 étant le "0" et signifiant aucune corrélation. Par ailleurs, un score inférieur à 1 implique une corrélation positive et un score supérieur à 1 une corrélation négative.
geary(x = data_immo$prix_med, 
      listw = neighbours_epci_w,
      n = length(neighbours_epci), 
      n1 = length(neighbours_epci)-1, 
      S0 = Szero(neighbours_epci_w))

$C
[1] 0.3324317

$K
[1] 18.24222

On voit dans les deux cas qu’il y aurait bien une auto-corrélation spatiale. Cela implique deux choses très importantes :

• La condition d’absence d’autocorrélation de nos résidus n’est pas vérifiée, le modèle classique n’est pas interprétable en l’état.
• La dimension spatiale joue un rôle, nous pouvons justifier d’étudier de manière plus approfondie l’autocorrélation spatiale et faire l’usage de la GWR.

data_immo$res_reg <- mod.lm$residuals

# Définition d'une palette de couleur
cols_v1 <- c("#08519c", "#3182bd", "#6baed6", "#9ecae1", "#c6dbef", "#eff3ff", "#ffffce", "#fee5d9", "#fcbba1", "#fc9272",  "#fb6a4a", "#de2d26")

# Réalisation de la carte
mf_map(x = data_immo, 
       var = "res_reg", 
       type = "choro", 
       border = "#ebebeb", 
       lwd = 0.1, 
       breaks=quantile(data_immo$res_reg,seq(0,1, by=1/11)), 
       pal=cols_v1, 
       leg_title = "Résidus de régression\nlinéaire 'classique'", 
       leg_val_rnd = 1)
mf_title("Résidus modèle lm") #titre

0.6.5 Petit bonus

Vous pouvez vérifier en une fois toutes les conditions de votre moèle de régression avec la fonction check_model()

check_model(mod.lm)